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Stabilization of polynomial dynamical systems 
using linear programming based on Bernstein polynomials 
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Abstract —In this paper, we deal with the problem of 
synthesizing static output feedback controllers for stabilizing 
polynomial systems. Our approach jointly synthesizes a Lya¬ 
punov function and a static output feedback controller that 
stabilizes the system over a given subset of the state-space. 
Specifically, our approach is simultaneously targeted towards 
two goals: (a) asymptotic Lyapunov stability of the system, 
and (b) invariance of a box containing the equilibrium. Our 
approach uses Bernstein polynomials to build a linear relaxation 
of polynomial optimization problems, and the use of a so-called 
“policy iteration” approach to deal with bilinear optimization 
problems. Our approach can be naturally extended to syn¬ 
thesizing hybrid feedback control laws through a combination 
of state-space decomposition and Bernstein polynomials. We 
demonstrate the effectiveness of our approach on a series of 
numerical benchmark examples. 

I. Introduction 

The problem of designing stabilizing controllers for non¬ 
linear dynamical systems is of great importance. In this 
paper, we study the problem of synthesizing static output 
feedback controllers for polynomial systems by solving 
a polynomial optimization problem to directly obtain the 
controller along with the associated Lyapunov functions that 
yields the proof of stability. 

Our approach inputs the description of a polynomial 
system and a desired region R to be stabilized. It then 
proceeds to find a static output feedback control law and an 
associated Lyapunov function to ensure local stability in R. 
Simultaneously, we ensure that the region R is an invariant 
of the resulting closed loop system. Our approach assumes a 
given structure for the feedback as a polynomial function 
of the outputs of the system. Furthermore, we assume a 
polynomial template form for the unknown Lyapunov func¬ 
tion. We proceed to encode the conditions for the Lyapunov 
function, obtaining a hard polynomial optimization problem 
that involves the coefficients of the Lyapunov functions and 
those of the feedback. 

The second part of the paper iteratively solves this op¬ 
timization problem through an iterative method variously 
called “V-K” iteration [19] or policy iteration [18]. The 
i th iteration of the approach selects a positive definite 
polynomial Vi and a feedback law Ui. Ideally, we require 
V( to be negative definite inside the region R for V to 
be a Lyapunov function guaranteeing asymptotic stability. 
Failing this, we first search for a new positive definite 
polynomial Vi +1 whose Lie derivative V( +1 has a larger 
maxima inside R fixing Ui , and adjust to a new feedback 
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law Ui +1 that improves the maximal value of V' +1 inside R. 
Each iteration is reduced to solving a Linear Programming 
(LP) problem using Bernstein polynomials combined with 
a reformulation linearization technique [5]. It is well-known 
that policy iteration does not necessarily converge to a global 
minimum, in general. However, our evaluation over a wide 
variety of benchmark examples shows that our approach is 
effective at converging to a global minimum by discovering 
an appropriate feedback law u* and an associated Lyapunov 
function V*. 

Automatic static output feedback design, or more gen¬ 
erally, finding feedback that satisfies given structural con¬ 
straints is well-known to be a hard problem in general. In 
fact, static output feedback stabilization of linear systems 
yields bilinear matrix inequalities (BMIs) rather than LMIs. 
A direct approach given by Henrion et al. [6] uses the 
characteristic polynomial of the transfer function matrix, 
and derives constraints that ensure the Hermite stability 
criterion for this matrix. As a result, they obtain a system of 
PMI (polynomial matrix inequalities), that is solved using a 
local optimization solver (PENBMI). In contrast, an indirect 
approach reduces the non convex BMIs to a series of convex 
LMIs. This was proposed as the so-called V — K iteration 
was proposed by El Ghaoui and Balakrishnan [19]. The ap¬ 
proach iteratively solves a bilinear problem by fixing one set 
of variables while modifying the other to result in a decrease 
in the objective values. The iteration alternates between the 
two sets of variables, until reaching a feasible solution. Our 
goal is to use this technique for polynomial systems while 
replacing BMI and LMI with linear and bilinear programs 
that can be solved more efficiently. A similar idea for solving 
bilinear problems appears in the work of Gaubert et al. [18], 
for finding invariants for discrete-time systems. Therein, the 
idea is called policy iteration. In this work, we will call 
our approach policy iteration , as well. The main differences 
between our work and that of El Gahoui et al. lie in our focus 
on polynomial systems, yielding more general polynomial 
optimization problems that involve the “V” variables relating 
to the Lyapunov function and the “K” variables relating to the 
feedback. Yet, by using policy iteration, we can separately 
focus on problems with a single set of variables at a time and 
use linear programming relaxations through a combination 
of Bernstein polynomials and reformulation linearization, 
discussed in our earlier work [5]. 

Existing approaches to stabilizing polynomial systems rely 
on linearization around the equilibrium. However, lineariza¬ 
tion can sometimes fail to be controllable, or yield region 
of stability that is much smaller than desired. Furthermore, 




Fig. 1. Overall structure of the controller synthesis problem considered. 

the output feedback stabilization for a linear system (or 
finding a feedback law satisfying a given structure) yields 
non-convex problems that are no easier to solve. Another 
class of methods (more related to our work) consists on 
reducing the problem to a set of LMIs or Sum-Of-Squares 
(SOS) formulations (see [20], [21] and references therein). 
In [21], an iterative SOS approach is proposed. This approach 
uses the Schur complement to produce a set of BMIs relaxed 
to an SOS problem. More precisely, an additional design 
nonlinear term e(x) is introduced, and causes bilinearity. 
An iterative approach is then obtained by fixing a guess 
for e(x) and iteratively updating it until feasibility is ob¬ 
tained. Once again, the major problem arises from the fact 
that the Lyapunov function and a static output feedback 
are needed simultaneously. Other approaches to controlling 
polynomial systems include the use of nonlinear optimal 
control techniques, feedback linearization, backstepping, and 
exact linearization. However, these techniques rely on the 
system being of a certain form and mostly involve state- 
feedback. A detailed comparison of the relative advantages 
of the direct approach presented here with other approaches 
to nonlinear stabilization will form an important part of our 
future work. 

II. Problems formulation and polynomial 

OPTIMIZATION PROBLEMS 
A. Problem formulation 

In this work, we consider a nonlinear control-affine system 
subject to input constraints : 

x{t) = f(x(t )) + g(x(t))u(y(t)), ueU. 
y(t) = h(x{t)). 

wherein x G M n represents the state variables, u G U 
represents the control inputs ranging over a compact set 
U CM P , and y G M 9 are the outputs. 

We assume that the functions / : R n -A M n , h : M n -A R q 
and the control matrix g : M n x M m -A R^ nxp "> defining 
the dynamics of the system are multivariate polynomial 
maps. The set of inputs U is a convex compact poly tope: 
U = {u G R p \ a u ,k ■ u < V& G K,u} where a u ,k G 
W, Pu,k £ M and Ku is a finite set of indices. Finally, we 
assume that x* = 0 n is an equilibrium for the system 0 . 
i.e /(0 n ) + g(0 n )u(h(0 n )) = 0 n . 

We define a region of interest R as a hyper-rectangle, R : 
[xi,xi\ x • • • x [x n ,x^\ with Xk <xf for all k G {1,..., n}. 

Stabilizing Feedback: In this work, we assume that the 
desired feedback is given by a function u : W -A U mapping 


outputs y to control inputs u to yield a closed-loop system 

X = f(x) + g(x)u(y), y = h(x) (2) 

We require that the closed loop system 0 be asymptoti- 
cally stable in R. This is achieved by ensuring two important 
properties. 

Problem 1 (Existence of Local Lyapunov Function): The 
system ([2]) has a local Lyapunov function V(x) in the region 
R such that 

1) V(x) is positive definite over R , i.e, V(x) >0 for all 
x G R \ {0 n } and V (0 n ) = 0. 

2) = ' (/( x ) + 9 { x ) u (h(x))) is negative definite 

over R. 

As such, a local Lyapunov function inside R guarantees that 
the system 0 is asymptotically stable in some neighborhood 
N of 0 n , where NCR. Specifically, N contains the largest 
sublevel set of V inside R as the stability region, but does 
not have to include R. To ensure that the system is stable 
inside all of R, we additionally require positive invariance 
of R. 

Problem 2 (Positive Invariance of R): The system ([2} is 
R-invariant , iff all trajectories with x(0) G R satisfy x(t) G 
R for all t > 0. 

Finding a feedback u(y) that solves problems [I] and [2] 
ensures asymptotic stability in the whole region R. 

Feedback Structure Finally, we consider feedback func¬ 
tions that conform to a given fixed structure. In other words, 
we consider feedback functions of the following form 

u(y) = H(y) • 0 = U(x) • 0 

where 6 G M 9 is a set of gain parameters to be determined 
by the synthesis procedure, the matrix H : M n -A 
is a given multivariate polynomial map that specifies the 
controller structure. Often, H is specified to include all 
monomial terms up to a given degree. However, more 
complex situations such as decentralized control may involve 
choosing specific structure for H. Figure [T] depicts the 
structure of the controller schematically. 

Let H(x) : H(h(x)) be the equivalent map as a function of 
the state variables. The input constraints (i.e. for all x G R, 
u G U) is then equivalent to 

Vfc G JCui Vx G R , au,k * R{x)0 < Pu,k- (3) 

Let O represent the values of 6 that satisfy Eq. 0. Under 
these assumptions, the dynamics of the controlled system 0 
can be rewritten under the form 

x(t) = f(x(t)) + G(x(t))8, 

where the matrix of polynomials G(x) = g(x)'H(x ), and 
<9gO. 

B. Reduction to polynomial optimization problems (POP) 

The first step is to fix a template form for the Lyapunov 
function V. We assume a polynomial form: 

v = v c (x) = c » x ° . 

\a\<D 











where a G N n , \a\ = JE a i> c : i c a)\a\<D are the unknown 
coefficients of the Lyapunov function and D e N is the 
maximal degree. 

We now focus on solving Problem [T] For a relatively small 
e > 0, this problem can be formulated as follows: 

1) Find a feasible set C s.t 

C : {c | minV c (x) — e||x|| 2 > 0} 

xGR 

2) Find feasible sets c G C' and 0 G O' s.t forall c G C’ 

and 0 G O', 

min —W c • (/ + GO) — e\\x\\ 2 > 0 

x£R 

Recall the set O from Eq. 0- 

Theorem 1: If C P| C' 0 and Of) O' / 0, then each 
c* G C P| C and 0 * G O P| O’ solves the local Lyapunov 
function existence problem (Problem [TJ. 

Proof: It is easy to see that the first condition will 
imply that V c will be positive definite, the second one implies 
that its derivatives ^ is negative definite. The last condition 
implies that the controller is admissible i.e u eU. ■ 

To solve the invariance problem (Problem [2]), we should find 
a controller (i.e a coefficient vector 0) ensuring that all the 
facets of the rectangle R are blocked. 

Definition 1 ( Blocked Facets ): A facet F of the hyper¬ 
rectangle R is said to be blocked for the system 0 if and 
only if 

\/x G F, np.{f{x) + G(x)0) < 0, 

where rip is its outer normal of the facet F. 

Let T denote the set of facets of the rectangle R , then solving 
Problem [2] can be formulated as follows : 

• Find feasible set Op such that for all 6 G Op s.t 

min np.(f(x) + G(x)0) < 0 , 

xeF 

for all facet F G T. 

Recall that O represents the feasible set from ([3]). 

Theorem 2: If Op P| 0/0, then each 0 * G Op P| O 
ensure the invariance of the rectangle R and solve Problem [2] 
where Op = [^j Op. 

FeR 

Proof: Since 0 * G Op then all the facets of R are 
blocked implying its invariance. The fact that G O' proves 
that the controller is admissible. ■ 

III. Reduction to Linear and bilinear feasibility 

PROBLEMS 

In this section, we are going to relax the previous poly¬ 
nomial optimization problems to a set of linear and bilinear 
feasibility problems. For doing so, we will briefly recall a 
relevant result showing how a general POP can be realxed 
to a linear program using Bernstein polynomials [5], then 
we will use this relaxation to build our linear and bilinear 
feasibility problems in order to solve our two given problems. 


A. Linear relaxation of a POP using Bernstein polynomials 

In this section, we are going to use Bernstein polynomials 
to establish lower bounds for our polynomial optimization 
problems (POP). More precisely, we seek tight lower bound 
for the optimal solution of the following POP: 

minimize p(x) s.t. x G R . (4) 

where p is multi-variate polynomial of degree S : 

(<5i,... ,S n ). 

We build a linear relaxation for problem Q, as follows: 

1) Change of variable qu mapping R to the unit box V = 
[0, l] ra . Let Pu = p O qu . 

2) Write pu in the Bernstein basis. 

3) Write an equivalent POP in the Bernstein basis. 

4) Exploit properties of Bernstein polynomials to formu¬ 
late a linear programming problem whose optimum is 
guaranteed to lower bound the POP in Eq. <0- 

We now explain the procedure in further detail. First of 
all, the mapping qu from any rectangle R to the unit box 
[0,1] n is an affine transformation. Therefore, the multi¬ 
variate polynomial pu is also of degree S and we can write: 

pu(y ) = for all y G U, 

a<5 

where (p a )a<s denotes the new coefficients of pu in the 
standard monomial basis, and the order relation a < 5 is 
such that OL{ < Si for all i G {1 ,,n}. By writing pu in 
the Bernstein basis we obtain the following form: 

Pu(y) = 

I<8 

where Bernstein coefficients (bi^)i<8 are given as follows: 


bps = 


E 

j<i 


n 

ji 


Si 

ji 


S n 

jn 


PJ = 


E 

j<i 


Pj- (5) 


and Bernstein polynomials are as follows: 


Bi,s(y) = 


y T (l n -y) 


5-1 


( 6 ) 


where y 1 = (yp,. ■ - ,y n ln ), 6 -1 = ( 61 -ii,... , 6 n -i n ) 


and 


Si 

k 


For the third step it is sufficient to replace the canonic form 
by the Bernstein form in the optimization problem, we then 
get the following optimization problem: 


minimize E bpsBps{y) 

I<5 

s.t y G U. 

zi = Bps(y). 


(7) 


The final step is now to remove the nonlinearities caused 
by the Bernstein polynomials by replacing each Bernstein 
polynomial Bj s by a fresh variable zj. In effect, we drop 




the relation zj = Bj^. To recover precision, we add some of 
the known linear relations between Bernstein polynomials: 

• Unit partition: = 1. 

I <5 

• Bounded polynomials: 0 < Bi,s(y) < 

B ItS {\), for all I <5. 

By injecting these properties in <(7j, we obtain the following 
linear relaxation: 

minimize E bl,6Zl,6 

I<6 

s.t z It s 6l, I <8, 

0 <z ItS <B ItS ($), I <5, [ > 

E z/ . 5 = 1 ’ 

I<5 


Lemma 1: The optimal value of 0 gives a lower bound 
for the POP 0. 


B. Linear and bilinear feasibility programs for existence of 
Lyapunov function (Problem^ 7]) 


Let V(x,c) be the assumed polynomial form for the 
Lyapunov function with unknowns c. We first focus on 
encoding the positive definiteness of V inside R. We recall 


the sets (7, C", O, O' from section [ILB 
First, we consider the set 



Let m{x) represent a vector of monomials involved in 

e||x|| 2 — V(x,c) so that we may write c\\x \\ 2 — V(x,c) : 

(fLm, where c = ^ ^ ^ for a suitable matrix L. Writing m 

in the Bernstein basis, we obtain m : Bz where z represents 
a vector of polynomials in the Bernstein basis and B is a 
linear transformation. Therefore, the problem ([8]) is written 
equivalently as 

min — c t LB z 

s.t. Az < b (9) 


Let C be the set of all values of c such that problem 0 with 
c £ C yields a non-positive optimal value. In other words, 

C : {c | (V z) Az < b => -c l LB z < 0} . (10) 

Lemma 2: C C C 

To represent the set C , we use Farkas lemma, a well known 
result in linear programming, to dualize eq. ( fT(j| ) and obtain 
our first linear feasibility problem for computing C C C. 

Lemma 3: The vector c is a solution to the problem 
in eq. © if and only if there exist multipliers c, A such 
that 

A f \ = -B t L t c, 6* A < 0, and A > 0 (11) 

Next, we consider the set O encoding the input constraints 
in 0- Let Hi denotes the Bernstein matrix associated to the 
i-th row of the the polynomial matrix H after mapping it to 
the unit box U (with respect to the degree S £ N n equal to 
the maximal degrees of H). Consider the set O defined as 


the feasible values of 0 that satisfy the following constraints 

<*u,k • BiO < /3u,k, & C ICu, Vi = 1,..., m . (12) 

Lemma 4: O C O. 

Now we will show that finding the feasible sets C’ and O' 
leads to a bilinear program. First, we can find a polynomial 
matrix B(x) to allow us to write 

-W c (*) ' (/(*) + G(*)0) = c* B{x)6, 

where 6 = ^ ^ ^ and B(x) = (W m (x))‘ • (/(x) G(x)). 

Here (VU m (x)) denotes the matrix where each column 
corresponds to the Jacobian of one of the monomials of the 
Lyapunov function. 

The main difference with the previous case is that instead 
of the vector of monomials m we have B{x)0. The degree 
S will be chosen as the maximal degrees of the polynomials 
in B(x). By consequence, the Bernstein conversion matrix 
will be a set of n matrices Be : i = Bi6 where Bi is the 
Bernstein conversion matrix corresponding to the polynomial 
row Bi(x) of the polynomial matrix B(x) after mapping it to 
the unit box U. Now using the same ideas as previously we 
will get by applying Farkas lemma a set of linear programs: 

Lemma 5: c is a solution to the problem in eq. ( fT0| ) if and 
only if there exist multipliers c and A 2 such that 

A t X L = —Bq/c , b f A 2 < 0, and A 2 > 0, for all i = 1,..., n. 

(!3) 

Since Bq^ = B[0, the previous lemma give us a set of 
bilinear feasibility problems for the feasible sets C’ and O'. 
But checking feasibility and solving a bilinear program is 
well-known to be NP-hard [17]. Rather than solve these 
problems directly, we consider a policy iteration approach 
in Section ITvl 

C. Linear feasibility programs for positive invariance (Prob¬ 
lem^ 

We now turn to the problem of encoding the invariance of 
the region R. Our approach reuses ideas from earlier work 
by Ben Sassi and Girard using the blossoming principle 
to enforce the invariance of a polytope for a polynomial 
system [4]. We obtain linear constraints over 0 that define a 
feasible region Of E Op such that choosing any 0 e Op 
guarantees that the region R will be maintained invariant. 

First, will need to define a facet and its outer normal [1] 
for a general rectangle R n = E&*A]: 

• 6c : {&/c>6c} | —^ {0,1} when for all k £ {l,...,n}, 
6c(a/c) = o and 6c(6c) = 1. 

• Fj£.( w .} = {x G R n | Xj = Wj }: the set of facets of 
R n where for all j E {1,..., n}, Wj £ {dj,bj}. 

• nj£.( w ) = (—l)(0(^)+ 1 )e J -: the outer normal of the 
facet Fj^y w .\ where the vectors ej form the canonical 
basis of R n . 

For the invariance context, all the results are derived from [4] 
so they are given without demonstration. We simply adapt 
the main result (Theorem 6 in [4]) to the specific form of 




the controller required in this work. For doing so we define 
for a fixed degree S = (5i,..., S n ), for all j G {1,..., n} 
and all / G {1,, Sj} : 

Ijj = {/ = (zi,..., i n ) G N n , such that / < S and ij = /}. 

More precisely, we need to replace in [4] the vector field / by 
f + G6 and the blossom values by the Bernstein coefficients. 
Let fu and Gjj denote the polynomial vector field / and the 
polynomial matrix G after mapping them to the unit box U 
and let fuj and Gup the associated Bernstein coefficient 
vector and matrix for all multi-indice I < 5. We will obtain 
the following result: 

Corollary 1: For all j G {1,..., n}, we have: 

1) The facet Fj^. ( aj ) of the rectangle R n is blocked for the 
controlled system x = / + GO if fu,i,j + GupjO > 0 
for all I G /j,o- 

2) The facet Fj^p b p of the rectangle R n is blocked for the 
controlled system x = / + GO if fu,i,j + Gu,i,jO < 0 
for all I G Ij : Sj , 

where fu,i,j and Gupj are respectively the j component 
(row) of the vector fup (matrix Gup)- 
The corollary gives us a linear program allowing to compute 
the feasible sets Of for all facets F G T. 

IV. Joint synthesis of polynomial Lyapunov 

FUNCTIONS AND CONTROLLERS 

First of all, we are going to present an algorithm to solve 
our stability problem, then we will show how the results can 
be improved by using a decomposition criterion and extend 
the results using this decomposition to a particular class of 
hybrid systems. 

A. Algorithm 

In this section, we give an algorithm allowing to summa¬ 
rize the previous results in order to solve our stabilization 
problems by synthesizing jointly the controller that stabilize 
the system and the Lyapunov function for the controlled 
system. 

In fact, the main problem when regrouping the feasibility 
problems of the previous section is that we have to deal 
with a bilinear program for which there is no practical way 
to solve it. We will define an iterative approach where for 
each step one of the parameters (0 for the controller or c for 
the Lyapunov function) is fixed and the other is computed 
by solving a linear program. The overall approach is given 
as follows: 

1) Initialize = 0. 

2) Compute feasible set C using feasibility problem GD- 

3) Find a ’’maximal” coefficient vector c G C for the 
Lyapunov function: 

We fix 0 = 0* and we solve the feasibility problems 
( p~3] ) by relaxing ” < 0” by ” < £” where t will be a 
positive decision variable to be minimized. The outputs 
of the linear program are (c*,£*). 

4) Find a ’’maximal” coefficient vector 0 for the controller: 
We fix c = c* and we solve the feasibility problems 
given by the (RHS) of (l2| ) and the ones of Corollary [I] 


By using the same idea of relaxing ” < 0” by ” < t” 
for a positive decision variable t and minimize over t , 
we get outputs (#*,£*). If t* ~ 0 STOP , else Go back 
to the previous step. 

When the algorithm terminates, the outputs (c*,0*) will 
give us the admissible controller and the Lyapunov function 
proving the asymptotic stability of the controlled system. 
The invariance problem of the rectangular domain will be 
ensured. 

B. Decomposition and generalization for a particular class 
of hybrid systems 

As mentioned in [5], the Bernstein relaxation ([8} can be 
much more efficient once a good decomposition is provided. 
By ’’good” we mean a box decomposition where local 
minima will belong to the edge of the box. Since the global 
minimum of the Lyapunov function is known in advance (0 n 
in our case), a decomposition of the rectangle R around zero 
(by putting zeros on the edges of the resulting rectangles) 
will significantly improve the precision of the approach. 
The drawback is that 2 n decomposition are needed. In fact 
by using this decomposition, each feasibility problem in 
the previous algorithm (except the invariance ones) will be 
replaced by 2 n feasibility problems. 

Now, since the approach deals with a box partition of the 
state space, one can easily extend the dynamical system 0 
to the following class of hybrid system where the state space 
is decomposed to boxes and each box has its own polynomial 
dynamic. More precisely, for all i < 2 n , let Ri be the set 
of boxes of our ’zero’ decomposition and the hybrid system 
will be following : 

{ ±i(t) = fi(x(t )) + G(x)0i , 0i e O X e Ri- (14) 

The difference here is that each of the 2 n feasibility problems 
(Step 4) will provide an admissible controller 0i trying to 
make the Lyapunov function decreasing in the corresponding 
box. So we will get a common Lyapunov function having 
multiple derivatives (one for each box). Also we should 
remark that when dealing with the invariance problem, linear 
feasibility problems of Corollary [I] should be adapted. In fact, 
for each box one should ensure the feasibility problems with 
respect to the facets that should be blocked. 

Remark 1: The previous result will hold for each other 
box decomposition. In fact we can always be reduced to the 
previous case by decomposing each sub box containing 0 n 
into sub boxes where 0 n will belong to the edges. 

V. Numerical results 
A. Illustrative example 

To illustrate the approach, we consider the following 2- 
dimensional polynomial system and a box R — [— 1, l] 2 . 

X\ — fl(x) = X 2 — x\ + Zx\ — 2 xi ^ 2 , 

= f 2 (x) = —Xi - Zx\ + x\ + 2 xix 2 . 

By simulation, one can see that the origin is not asymptoti¬ 
cally stable and that the box [—1, l] 2 is not invariant for the 




Fig. 2. Vector fields and some trajectories of the uncontrolled system. 

Fig. 3. Vector fields and some trajectories of the controlled system. 


system (see Figure [3]). Using our approach, we aim to find 
a linear state feedback controller ensuring the asymptotic 
stability of the origin and the invariance of R. We will 
consider the following controlled system: 


x(t) = f(x(t)) + g(x(t))u(x(t)), 


where g{pc) = I 2 = 
A = 


and u{x) = Ax where 


1 0 
x 0 1 

an <212 

a2i a22 

Since we look for a linear state feedback controller, we can 
write u{pc) = H(x )6 where 

X\ X 2 0 0 

0 0 X\ X 2 

(an, ai2, <221, <222) • 

For the Lyapunov function, we fix the following form : 


H(x) 


and 0 


V c (x) = C\X\ + c 2 x 2 + c 3 x\ + c±x\ + C 5 X 1 X 2 + c 6 xf + cix\. 


We impose that — 5 < q < 5 for all i G {1,..., 7} and add 
the fact that c* > 0.01 for all i G {3,4} in order to ensure 
that V is positive definite. Also the linear coefficients of the 
controller are bounded by —5 and 5. 

The iterative approach needs two iterations to globally sta¬ 
bilize R. Outputs are : 

f -4.5471 0.7000 \ 

9 A ~ \ 3.9290 -4.6218 )' 

. V(x) = 0.01(x?+^)+0.009xix 2 +0.036xf+0.023x|. 
One can simulate the obtained system and verify the asymp¬ 
totic stability and the invariance of R (see Figure [3]). Now, 
we will use the approach to deal with the Hybrid case. More 
precisely, we decompose R around zero (Ri = [—1,0] 2 , 
R 2 = [-1,0] x [0,1], R 3 = [0,1] x [-1,0], R 4 = [0, l] 2 ) 
and try to find for each sub-box Ri a linear controller ui 
such that the following Hybrid system 


x(t) = f(x(t)) + g(x(t))ui(x), 

is globally stable with respect to R where Ui(x) = Aix for 
all x G Ri and all i G {1,..., 4} . 

In this case, only one iteration is needed to stabilize the 
system inside R since we have more freedom in the choice 
of the controller. Outputs are : 


-4.4721 

-3.4219 

-2.9376 

-4.0957 

-4.3795 

0.3130 

1.1904 

-4.3770 

-4.3331 

2.6016 

3.3924 

-4.2926 

-4.1427 

-3.0418 

-3.3052 

-4.4195 


• A\ - 

• A 2 = 

• A 3 = 

• A4 = 

. V(x) = 4.7737xf + 4.7743x1 + 4.8172xf + 4.8175x 


By simulating trajectories in those boxes, we can verify 
that the stability property and the box invariance hold (see 


Figure]?] for R 1 and R 2 ). 



Fig. 4. Vector fields and some trajectories of the controlled system 
associated to R\ (on the top) and R 2 (on the bottom) 











































B. Benchmarks 

We discuss the results obtained for a set of benchmarks 
borrowed from the literature. We run the algorithm until a 
good precision e Qs reached or a fixed number of iterations 
(the approach fails to make progress). In the latter case one 
can add more flexibility in the templates by adding terms of 
higher degrees. In failure cases, we remove the invariance 
constraints in order to achieve just the asymptotic stability 
property. We report separately stability (Stab column) and 
invariance (Inv column). A threshold of precision around 
10 -6 is considered to confirm that the property holds. 
We report also the number of iteration needed to achieve 
the given precision. A detailed description of the systems, 
explicit expression of Lyapunov functions and controllers are 
given in the Appendix. 


TABLE I 

Table showing performance of our method on a set of 

BENCHMARKS. 


Id 

R 

u 

e 

Stab 

Inv 

Iter 

1 

[— 0 . 5 , 0 . 5 ]^ 

[-1,1] 

4 * 10 “ ly 

/ 

X 

1 

2 

[-M] 2 

[- 2 , 2 ] 

4 * 10“ 9 

/ 

/ 

2 

3 

[- 1 , l ] 2 

[- 4 , 4 ] 

2 * 10- 17 

/ 

X 

1 

4 

[-M] 2 

[-1,1] 

4 * 10- 7 

/ 

/ 

3 

5 

M ,1] 3 

[- 10 , 10 ] 

2 * 10“ 6 

/ 

/ 

4 

6 

[— 0 . 5 , 0 . 5] 3 

[- 5 , 5 ] 

2 * 10- 7 

/ 

X 

6 

7 

[— 0 . 5 , 0 . 5] 3 

[-3,3] 

9 * 10- 7 

/ 

X 

3 

8 

[— 0 . 5 , 0 . 5] 3 

[-1,1] 

4 * 10“ 5 

? 

X 

9 

9 

[— 0 . 1 , 0 . 1] 4 

[- 5 , 5 ] 

5 * 10“ 5 

? 

X 

3 

10 

[— 0 . 1 , 0 . 1] 4 

[- 10 , 10 ] 

8 * 10“ 5 

? 

X 

4 

11 

[— 0 . 05 , 0 . 05] 5 

[-1,1] 

6 * 10“ 5 

? 

X 

2 


Note that that invariance conditions usually make the 
feasibility of the approach very restricted since it needs 
to holds simultaneously with the stability conditions. This 
explains the fact that only few stabilazable systems can only 
have the invariance box property. The computation time is 
roughly in size of the problem and the templates: roughly 
each iteration of two dimensional systems (systems 1, 2, 3,4) 
required almost one second, for three dimensional systems it 
required between two and three seconds (systems 4, 5, 6, 7). 

VI. Appendix 

Example 1: (see [22]) 

f x = y. 

\ y = —x + u(y). 

• u(y) = —2 y. 

• V(x,y) = 0.01(x l 2 + y 2 ) 

Example 2: (see Lectures on back-stepping 

J x = y — x 3 . 

\ y = u{x,y). 

• u(x, y) = -x - \y + |a; 3 . 

• V (x, y) = 0.01(y 2 + x 2 y 2 ) + 0.0102a; 2 + 0.0007a;2/. 


l e denotes the precision t* of the algorithm. 

2 h ;tp://control.ee.ethz.ch/ apnoco/Lectures2014 


Example 3: 

f x = y 

\ y = u(y)y 2 - x. 

• u(y) = 4(t/ 2 - y). 

• V(x, y) = 0.01(x 2 + y 2 + x 2 y 2 ) + 0.005(x 4 + y 4 ). 
Example 4: (See [23]) 

{ x — —x(0.1 + (x + y) 2 ) 

\ y = (u(x) + x)(0.1 + (x + y) 2 )- 

• u(x) = —X. 

• V(x, y) = 0.01 (y 2 + x 2 y 2 ) + 0.0657x 2 + 0.0022 xy + 

O.OOlV. 

Example 5: 

{ x = y + 0.5 z 2 . 
y ~ z. 

z = u(x,y,z). 

. u(x,y,z) = -0.59185x - 5.9217 y - 0.51825 2 ? + 
0.061785x 2 + 0.12415 xy - 0.4642^ + 0.048453x 3 - 
0.57345T/ 3 . 

. V(x,y,z) = O.Olx 2 + 0.0583 t/ 2 + 0.0099z 2 
0.0134x7/ + 0.003x2: + 0.004 y 4 + 0.0024 yz + 0.00032: 
Example 6: (See [24]) 


{ x = —x + y — z. 

V = -x{z + 1) - y. 
z = —x + u(x, z). 

. u(x, z) = 1.76524X - 4.7037^. 

• V(x, y, z) = 0.01(x 2 + y 2 ) + 0.0132: 2 . 
Example 7: (see Lectures on back-stepping) 

{ x = —x 3 + y. 

y = y 3 + z- 

z = u(x,y,z). 


• u(x,y,z) = —0.083339x — 3.5413 y — 0.338682: — 
0.4325x 3 . 

• V(x, y , z) = 0.01(x 2 + z 2 ) + 0.0333^ 2 + 0.0033x7/ + 
0.0048x2: + O.OO6I7/2:. 

Example 8: (See [24]) 


x = z 3 — y. 
y = z. 

Z = 7i(x, 7/, Z). 


• tx(x, 7 /, z) = —0.86597x — 0.162087/ — 0.615972:. 

• V(x, 7 /, z) = 0.01(x 2 + z 2 ) + 0.0333^ 2 + 0.0179x7/ + 
0.0129x2: + 0.0127 t/ 2 . 

Example 9: 

x = y- 

y = — 0.1t/ — IO2: + xv 2 . 
z = v. 

v = —2: — v + 7 i(x, 7/, 2:, v). 

. u(x,y,z,v) = -12.0271x - 8.1243t/ - 10.2755^ - 
10.047t;. 

. V(x,y,z,v) = 0.1202x 2 + 0.01(t/ 2 +t; 2 )+0.2201^ 2 + 
0.2556x2: + O.OIOIxt; + 0.01578t/2: + 0.0115t/t;. 


+ ^ ’ 












Example 10: (Ball and Beam example [25]) 

x = y- 

y = — 9.8z + 1.6z 3 + xv 2 . 

. 

Z = V. 

V = u{x). 

• u(x) = — 6x. 

• V(x,y, z,v) = 0.0672a; 2 + 0.01 y 2 + 0.1074z 2 + 
0.0136?; 2 — 0.0043x7/+0.149a; z+0.0023a;?;+0.008?/z+ 
0.0189 yv — 0.003z?;. 

Example 11: 

x = —O.lx 2 — 0.4a;?; — x + y + 3z + 0.5?;. 
y = y 2 — 0.5 yw + x + z. 

< i = 0.5z 2 + x — y + 2z + 0.1?; — 0.5??;. 
v = y 2z + 0.1?; — 0.2??; + u(x, y , 2 , ?;, ??;). 
w = 2 — 0.1?; + ??(x, y, z, v, w). 

• u(x, y , z, ?;, ??;) = — 1.5a; — 1.5?/ — 1.5 2 — 1.5?; — 1.5??;. 

• V(x,y,z,w,v) = 0.01(x 2 +y 2 +v 2 +w 2 ) — 0.0066xy — 
0.0252 xz — 0.008(o;?; + yv) + 0.005a;??; + 0.001 yz + 
0.0167 yw — 0.0023^?; — 0.0121^??; + 0.001?;??;. 

VII. Conclusion 

In this paper a linear programming approach is presented 
allowing to deal with the stabilization problem of polynomial 
systems. The approach is based on Bernstein polynomials 
and propose a policy iteration technique allowing to avoid 
bilinear programs by having an iterative approach of linear 
programs instead. The benchmarks results show that the 
method can be efficient in practice. The drawback of this 
technique is that no convergence result is guaranteed and 
even in case of convergence there is no guaranty that it will 
be to a local minima. A future work will be a deeper study 
of the failure case or the fix point (once the algorithm result 
does not improve): an idea is to fix small variation for each 
variable of the bilinear program and try to find a descent 
direction helping the algorithm to improve. 
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